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INTRODUCTION 


This Parameter Estimation Supplement describes the PEST computer pro- 
gram and gives instructions for its use in determination of lunar gravitation field 
coefficients. 

PEST was developed under Contract NAS 5-11900 for use in the RAE-B 
lunar orbiting mission as a means of lunar field recovery. The observations pro- 
cessed by PEST are short-arc osculating orbital elements. These observations 
are the end product of an orbit determination process obtained with another program. 

PEST's end product is a set of harmonic coefficients, C and S , to be used in 

nm nm 

long-term prediction of the lunar orbit. PEST employs some novel techniques in its 
estimation process, notably a "square-root” batch estimator and linear variational 
equations in the orbital elements (both osculating and mean) for measurement sen- 
sitivities. 

Section 1 describes the program's capabilities, giving operating instructions 
and input/output examples. This section supplements the User's Manual for the 
MAESTRO program (Ref. 1). PEST utilizes MAESTRO routines for its trajectory 
propagation. Section 2 describes PEST's program structure and contains descrip- 
tions of those of its subroutines which are not common to MAESTRO. The MAESTRO 
subroutines are described in the Programmer's Manual for MAESTRO (Ref. 2). 
Section 3 contains some of the theoretical background information for the estimation 
process. Appendix A contains a derivation of linear variational equations for the 
Method 7 elements. 
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I. PEST PROGRAM DESCRIPTION (User's Guide) 

1. 1 Highlights and Capabilities 

The parameter estimation (PEST) program is a weighted least squares 
batch processor which estimates 

a. coefficients of the lunar gravitational field 

b. initial orbital elements 

c. solar pressure coefficients. 

The measurements used in the estimation are short-arc determinations of the 
orbital elements relative to the Moon. PEST can process all or any subset of 
the residuals in the elements. The elements can be osculating or averaged 
values of either classical form or Gaussian. 

Methods 4 or 5: p , e , M , to , i , £2 

Method 7 : p , e sin m , e cos to , to + f , i , Q 

Method 8: p, e sin to, e cos to, tO + M , i , £2 

The measurement data can be 

a. read in from magnetic tape 

b. read in on cards 

c. simulated by the program itself 

Measurement sensitivities can be generated either by the secant method or by inte- 
gration of linear variational equations. The sensitivities may be "frozen" at some 
iteration of the estimation process, thus saving the computation time needed to re- 
compute sensitivities. The estimated state increment may be included with a gain 
factor other than unity to combat nonlinearity problems. Selected parameters may 
be held fixed during the estimation while they are being "considered" in the calcu- 
lation of the state increment. 
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1. 2 Operating Instructions 


Figure 1. 1 illustrates the order in which PEST accepts and expects input 
(card) data. The READ statements occur in only two subroutines, INPUTF and 
FILTER. INPUTF reads data into INPUT common, reading three (address, value) 
pairs per card with a (IX, 3(15, D15. 8)) format. A blank or zero in the first ad- 
dress field causes a return from INPUTF. 


The first or only card read by FILTER is to identify the parameters to be 
estimated and to flag the primary options of this capability. The next card is a 
measurement selector. FILTER next expects measurement data. If an end of file 
is encountered at this point, the (only) case will be executed. If another case is to 
follow the current one, an F10.2 zero must be read when FILTER expects a 
measurement. Control then is returned to the main program and subsequently to 
INPUTF before returning to FILTER. 


1 . 

2 . 

3. 

* 4. 

* . 5 - 
* 6 . 

** 7 # 

** 8 . 
** 9. 


INPUT common data 

Blank or 16 zero to exit from INPUTF 

IX array card(s) to identify state variables and primary options 
Measurement selector card 
Considered parameters card 
Measured elements and their epochs 

Blank or F10.2 zero to signal "data for next case following" 
INPUT common data for case 2 
Steps 2-7 for case 2, etc. 


* optional 

** needed only if more than one case is to be run 
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PEST Main 



Figure 1.1 


Input Data Read Logic 







INPUT Common Data to Consider 


1044 

1030 

1033 

1049 

1036 

4 

180 

1068 

1069 

450 

460 


1035 

1094 

1095 
1084 


50-55 


56-76 


MODE 

KOUT 

KRW 

KSAD 


KMETH 

TFINAL 

DELT 

NINT 

NORD 

TOUT 

STEP 


MODLEM 

NMOD 

MMOD 

KSOLP 


General 

Flag to invoke PEST. Set to 8 . 

Output flag. Set to -1 . 

Read-write flag for logical 22. Usually 1. 

Shadow flag. Set to 0 . 

Trajectory Propagation 

Propagation method flag. May be 4 , 5 , 7 , or 8 
unless KVAR = 1. 

Final time (not used in estimation, where 
measurements control stop). 

Integration step size (see KMETH). 

Number of intervals used in averaging. 

Number of ordinates used in averaging. 

Print table (set it larger than last measurement time) 

Time interval for simulated measurements or when 
KRW is zero. 

Field Model 

Lunar field model flag (3 for JPL 15, 8;1 for LI). 

Order of zonal harmonic model. 

Order of tes serai harmonic model. 

Solar pressure flag (must be nonzero if SP is estimated). 
Epoch 

Input epoch to which all times are referred. 

Measurement Weights 

Measurement weighting matrix (upper half of 
symmetric matrix). 
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Initial Conditions 


1019 

KINPT 

Input coordinate flag (see KAVST). 

30-35 

EL MI 

Initial elements (see KAVST). 

1098 

KAVST 

Averaging and initial condition flag. 

0 neither input elements nor measurements are averaged. 

1 average input elements if KMETH is 5 or 8 . 

2 average measurements but not initial elements. 


3 average both initial elements and measurements. 
-1 average nothing, but set initial elements to first 
measurement. 

-2 average measurements, set initial elements to 
first measurement. 


1096 

1064 

1086 


KVAR 

NGROPT 

KZIP 


Sensitivity Calculation 


Sensitivity method flag. 0 for secant, 1 for variational. 

Number of trials for which sensitivities are to be 
recomputed. 

Flag for ignoring implicit term in variational equations. 


reification of State Variables and Options 


The IX array is read in by means of one or more records (cards) written 
in a 2013 format. Reading of these records is terminated when the (N+l)st 
element of IX is negative or zero, or when five records (100 elements) have been 
read. Valid identification numbers are: 


Field coefficients IX( i) = 2-271 except MOD(IX( i),16) = 0 or when 

IX( i) lies outside the field specified by 
NMOD, MMOD. 

Initial conditions IX(i) = 273-278. 

Solar pressure IX( i) - 279, 280. 

The first nonpositive identifier controls the operational logic: 

IX(1) = 0 Generate simulated measurement data. 

IX(N+1) = 0 Estimate the state if N is greater than zero. 

IX(N+1) =-l Compute measurement residuals and (if N is greater than 0) 
compute/print measurement sensitivities without actually 
estimating a new state. 
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The field identifiers can be obtained with the assistance of Table 1.1. Find the ap- 
propriate coefficient in the table and add its row-number (shown on the right) to the 
number shown above its column. This sum is the identifier to use in the DC array. 
The order in which you specify these numbers is now immaterial — IX will be 
sorted appropriately by the program. 

Examples: 

0 Simulated data will be generated from time zero to TFINAL by STEP. 

2_-l Residuals and sensitivity to C will be printed. 

_34274 0 Estimate C ^ and either e Q or e Q cos , depending on KMETH. 

Measurement Selector Card 

The measurement selector (IM) card is read by means of a 2013 format. 

The first six words are the numbers (1-6) of the elements to be used in processing. 

A zero in the field will exclude the corresponding element. (1 , 2 , 3 , 0 , 5 , 6 excludes 
U-'+f residuals if Method 7, W+M if Method 8, is used, or w if Methods 4 or 5 
are used in propagation.) The next number is the number of observations to be pro- 
cessed. If this is greater than the number of observations to be found on tape, only 
the smaller number will be processed, of course. The next number dictates IQ , the 
weighting matrix treatment. If IQ is 1 , Q is I; if IQ is 2 , Q is constant; 
and if IQ is 3 , Q is variable and is read in from the tape with the observation. 

In this latter case, the input weighting matrix is assumed to describe the uncertainty 
in the Cartesian state and will be converted appropriately by the program to an ele- 
ment weighting matrix. If Q is to be treated as a constant, it must be read in 
(locations 56 through 76) as an element weighting matrix, as it will not be converted 
by the program. In either case, Q must be positive definite. 

ITMAX, the next indicator on the card, limits the number of iterations allowed 
in estimation if convergence is not achieved. The estimator will quit after passing 
through the data span ITMAX +1 times. The next indicator on the card is IGAIN, 
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0 

1(> 

32 

48 

G4 - 

80 

9G 

°i,o 

C I,1 

S l,l 

S 2,l 

S 3.1 

S, 1 , I 

S 3.1 

°2,0 

°2,1 

C 2,2 

S 2,2 

3,2 

*4,2 

K 5.2 

°3,0 

c :*,i 

C 3,2 

C 3,3 

S 3,3 

S 4,3 

S 5, 3 


°4,0 

*4,1 

°4.2 

°4,3 

C 4,4 

S 4,4 

s r>. 4 

C 5, 0 

C 5, l 

C 5,2 

e 5,3 

C 5,4 

C 5. 5 

8 s, 5 

C (i,0 

°G,1 

°G,2 

C G , 3 

°G , 4 

C G,G 

. C G,G 


c io,o 

c it,o 

C 12,0 

C 13,0 

C 14,0 


f 


Table 1. 1 


Field Identifier Aid 


112 

.128 

144 

ICO 

176 

192 

208 

224 

240 

256 


*6,1 

*7,1 

S 8,l 

S 9,l 

S 10,l 

s n,i 

S 12,l 

S 13, 1 

S 14,l 

S 15,l 

1 

*0,2 

*7.2 

*8,2 

*9, 2 

S 10,2 

*11,2 

8 12, 2 

S 13,2 

S 14, 2 

S 15,2 

2 

S G,3 

S 7, 3 

S 8,3 

S 9,3 

S 10, 3 

S ll,3 

S 12, 3 

S 13,3 

8 14,3 

S 15,3 

3 

S G, 1 

S 7,4 

S 8,4 

®9,4 

S 10,4 

*11.4 

S 12,4 

*13,4 

S 14,4 

S 15,4 

4 

s . _ 

S 7, 5 

*8,5 

S 9,5 

S 10,5 

S ll,5 

S 12,5 

S 13,5 

S 14,5 

S 15,5 

5 

s 0>, G 

*7,G 

S 8,6 

S 9,G 

S 10,6 

S ll,6 

S 12,G 

S 13,G 

fl4,6 

S 15,G 

6 

C 7,7 

*7,7 









7 


C 8,8 

S 8,8 








8 



°9,9 

S 9,9 
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C 10,10 

8 10, 10 
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C ll,ll 

s n.ii 
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C 12,12 

*12,12 




12 







C 13, 13 

S 13,13 



13 








C 14, 14 

S 14,14 


14 







i 

t 

C 15, 15 

*15,15 

15 



the reciprocal gain factor on the incremental state to be added to the state between 
trials. Thus, if IGAIN = 2 , only half of the estimated increment will be added. 
LTOL , the next number on the card,’ controls the tolerance on the iteration as fol- 
lows: TOL = 1. D-LTOL . The next and last number is an indicator for the 
''consider’' option. If ICON = 0 , nothing happens in this regard. If ICON ^ 0, 
another card is expected to identify the "considered" parameters. The "considered" 
parameters are identified by the same code as the IX array. 


Input of Measured Elements 


The current means for input of measured elements is by means of DATA5. 
Each record contains 8 words written in an 8F10.2 format. These words are: 


YYMMDD Year, month, and day of epoch (710824 is 8/24/71) 

HHMMSS Hour, minute, second UMT (183025.5 is 18:30:25.5) 

A Semi -major axis (km) 

E Eccentricity 

F True Anomaly (deg) 

W Argument of pericenter (deg) 

I Inclination (deg) 

O Longitude of the ascending node (deg) 


If KMETH is 5 or 8 and KAVST is -2, 2, or 3, the input measurements are • 
considered to be osculating and selenographic and will be averaged and converted to 
mean of 1950.0 before being written onto logical 22. If KMETH is 5 or 8 but KAVST 
is -1, 0, or 1, the input measurements (and epochs) will be considered to be al- 
ready averaged and converted to mean of 1950. 0, so will be simply converted to 
KMETH variables before being written onto 22. If KMETH is 4 or 7 (KAVST 
should then be -1, 0, or 1), the input measurements will be considered to be os- 
culating and selenographic and will be converted to KMETH variables, mean of 1950.0. 
The longitude parameter for KMETH = 8 is treated as a) + f outside the propagation 
routines when computing residuals and sensitivities. 
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Input Measurements from. GTDS 


The GTDS measurements are expected on magnetic tape in the form of 

YYMMDD Year, month, and day of epoch 

HHMMSS Hour, minute, and second of epoch 

Cartesian position components (km), 
selenocentric, Equator of 1950. 0 

X -| 

* Cartesian velocity components (km/sec), 

selenocentric, Equator of 1950. 0 

Z -* 

Definition of Peripheral Units 

The parameter estimation program uses logical units 21, 22, and 23 as se- 
quential datasets which must be defined in the JCL of the job setup. Logical 21 is 
used for input measurement data from GTDS. (Two words of epoch and six com- 
ponents of Cartesian state, EE50 in an 8D23. 15 format, requiring a record length 
of 188 bytes.) Logical 22 is used for storage of the measured time and elements. 

It therefore requires a record length (LRECL) of 60 (8*7+4) or greater and a block- 
size of J*LRECL+4. Logical 23 is used as a scratch dataset for storage of the 
weighted sensitivities. It requires a record length of 48*N+4 or greater, where 
N is the number of state variables included in the case. 



1. 3 Input/ Output Examples 

Figure 1. 2 shows the input which resulted in the output of Figure 1.3. 
The INPUTF data is abbreviated in Figure 1.2, but fully printed in Figure 1.3. 
The important numbers of the INPUTF data are: 
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FORTRAN Coding Form 



program 

PUNCHING 

INSlRUChONS 

GRAPHIC 





pkoghammip 

DATE 

PUNCH 






H 

MAlTMf Ml 
NUMpf P 

6 

u 

FORTRAN STATEMENT 

nr 

2 3 4^ 

6 



' ! 3 

0 

±23 

mm 

1 

jjjj 

H 


i: \-'ao 


ill! 

i H, 

I ! 'l 

32 7 23 

Xxt\ 

HHI 


IHS 



, 1 1 

) i ; 


1 

i 1 

i 

! ! 

! ( ' I 

1 I , i 

1 1 

: 1 j ! 

1 1 H : 

Mil 

! 1 U 

1 i ! | 

• 1 1 
ill; 

, j ! ; 

— 


wl 

1 

.or } 

HM 

S3 

Tj 

warn* 

' ~o/n 

-?f 1 1 i 

III! 
) 1 1 1 

i ' j 

1 


ill 1 

BHhHH 

i ! 1 i 

BHi 


i 

H 

1 

! ! 1 

1 

■ 

■ 

1 

— 

■ 

i 1 : 

1 1 i 

! II 1 

HI 

1 

1 


! ' 



i 1 

' J 

■HI 

Ml! 

HHI 


'■lO c ) 

% 

yu 

HI 

■ 

1 

HI 

M 

MHb 

4 '/*! 

H 

1 

1 

l 

i 

1 1 



i i 1 

ill! 

MM 

IH 

H 





i 

Q 

$ 

nr 


IH 


u 


ill 


1 

1 

1 

! I 
1 1 ! 

r | 

1 1 

MM 

HHI 

■HHHH 


X 7 

WEE 

3 

4 57) 

s 

H 

? 7mt 

35LJ&J2 

■HH 

7-20-3,2 

0 1:! 

7221 

727? \ 

B— 


i i ; 

1 : ; 


'/■ : 

% 

mm 

* i T 

■ ! 4 / 

a i ! h 

3 i !/ 

HHI 

j i i ! 

! ! i ! 

HHI 

i : i 


RfljE 

■ 

1 


in 

m 

HI 

H 

ran 

inn 

n 

HHI 

HHI 

IHHI 

HHHHH 

■ 

HH 


n 

■ 

1 

SB 

HI 

in 

HI 

n 

IHH 
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HH 

1 

■ 

in 

■ 

an 

IUHI 

■HHHH 

in 
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Figure 1.2 Estimation Input 




jriuuKH-. i. 3 


DATA CAR 

DS F3S CASE 

i 





- 

“ " 

30 

0.283 80 OC 000 

04 

31 

C.395700CCCD- 

02 

32 

0. I 2352 1 CO 00 

03 

33 

0. 1 C4699CC00 

03 

34 

C • 1 16156 COCO 

03 

35 

0*36399000 00 

02 

50 

0*6000000000 

01 

5 1 

C.3CCC00CCCD 

02 

52 

0. 1 973 00 00 OD 

04 

53 

0.1 5C COOOOOO 

02 

£4 

C. 14Q0C0CC0D 

02 

55 

C*43G0000QCD 

02 

1044 

o.accoooocoo 

01 

1030 

-C. lOCCOOOCOO 

01 

1 060 

0*0 


4 

0*3456000000 

06 

460 

C.43200000 OD 

05 

450 

0. 1 OCCOOOOOD 

08 

1035 

0.3C0C000C0O 

01 

10 9 4 

C* 40G000CC0O 

01 

1 095 

C. 4000000000 

01 

1398 

0*0 


0 

0*0 

. . 

0 

0*0 


1036 

0.8C0CCCC00D 

01 

1 60 

0.4320000000 

05 

1 C 98 

-0* 200000000D 

01 

1036 

0.7COCOOCOOO 

01 

1 60 

0 . 12 C C 00 0000 

04 

1 C 98 

•0* 1 OCOOOCO OD 

01 

10 96 
1096 

0*0 

0*1000000000 

01 

0 

1066 

C.O 

C. lOCCOOOCOO 

or 

0 

0 

0*0 

0*0 ~ ■ 


IN1T IAL 

JULI AN DATE 

2441664. 

135 2 

— 

— 

. 


— 

- 

INITIAL 

CONDITIONS 

— 


— - ' • 



IN EARTH MEAN EQUATOR AND EQUINOX OF 1950 
6/30/73 15 NRS 14 MINS A3. 000 SECS 


X Y Z 

0*745050720 03 0.2S7-J5274D 04 ' -0 . 9436 94220 03 

XOOT Y COT ZODT 

0*945352200 00 - C . £4 1 £6 C 760-0 1 -0.9C7426940 00 

ATTITUOE INPUT KEY ILCC 1C87) - 0 

CURRENT ATTITUDE. RA = -2.201 DEC = -43.780 ' 



IM 



M IQ IT IG LT 


12 

3 

4 

5 

6 10 1 318 



2 


C ( 

2* 

0) -0.1996C-C3 

2003 

1 8 


C { 

2 * 

1) 0*8171 C— 05 

20 19 “ “ 

49 


5 i 

2 t 

1 ) -0.72 13C-05 

23 19 

34 


C ( 

2, 

2) 0.2359C-C4 

2035 

- 50 


S ( 

2 « 

2) 0.4538C-05 

2335 

3 


C ( 

3 9 

0> — 0 * 5878C— 05 

2004 

1 9 


C ( 

3* 

1) 0.3CC1C-C4 

2020 

65 


S ( 

3. 

1 ) 0 . 1421 C— C5 

23 20 

35 


C ( 

3* 

2) 0.4698D-05 

2036 

66 

- • -- 

5 C 

3* 

2) 0*57480-06 

2336 

51 


C ( 

3* 

2) 0.4647C-C5 

20 5 2 

67 


s< 

3, 

3) —0*2919 C— C5 

2352 

273 




" P 


276 


— 


W+M ■ -~u/rf ~~ ' ~ 

- — — 

277 




INC 


273 

— 

— 


LAN - - - 




PARAMETER ESTIMATION 


TRIAL NO. 0 
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2072.151 
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730 6 14 • 

225641* 

-2047. 9C9 

1 890 .460 

-741 .893 - 

0*87024 -0*64524 

0*74686 

73061b * 

105641* 

-2298*148 

1670 »8C8 

-502.309 - 

0*73142 -0.75892 

0*78677 

7306 15 * 

22564 1* 

-2503*002 

14 17.462 

-252.365 - 

0*57812 -0*85729 

0*81129 

7306 16 • 

105641* 

-2658.356 

1135.349 

- 3.288 - 

0.41343 -0.93837 

0*81989 

730 6 1 6 * 

225641. 

-2760*966 

830 ,001 

259.812 - 

0*24071 -1*00049 

0. 81 238 

7306 17 i 

105641* 

-28C8.672 

507 .483 

512.230 - 

0.C6356 -1*04244 

0*78886 

7306 17 • 

22564 1. 

-2800*316 

174.295 

755.508 

0*11430 -1*06343 

0*74979 

7306 18 ♦ 

105641* 

-2736*05 1 

-162.767 

964.684 

0 * 28 9 1 C -1*06317 

0.69601 

730613* 

22564 1 * 

-261 7.298 

— 496 * 7? 7 

1 1 95.024 
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Figure 1. 3 (cont. ) 


1044 MODE is 8 to signal parameter estimation 

1030 Output key is negative to suppress extra output 

1060 Zero signals no print of input common 

1035 3 signifies the JPL 15,8 field is estimated 

1094 4 means the field is truncated at n = 4 

1095 4 means the field is truncated at m=4 

1036 Method 7 propagation is used (Gauss' elements, no averaging) 

180 The integration time step is 20 minutes 

1098 Initial conditions (time and elements) are taken from the 

measurement tape when KAVST = -1 

1096 KVAR is 1 to use linear variational equations for measure- 

ment sensitivities 

1086 KZIP is 1 , but has no significance when 1036 is 7 

The zero in column 6 of Figure 1. 2 causes a return from INPUTF and continuation 

to FILTER where the next card is read. This card is to select parameters to be 

estimated. For the case shown, these were all of the gravitational coefficients of 

a 3,3 field plus four initial elements. The significance of the numbers is indicated 

in the printout shown in Figure 1. 3. The next card says that residuals in all six 

Method 7 elements are to be processed in the estimation. It also says that only the 

first 10 sets of elements on the measurement file are to be processed, that no 

weighting is to be done (IQ = 1), that no more than three iterations are to take place, 

that unit gain is to be put on the included state increment, that the convergence 
-8 

tolerance is 10 , and that no parameters are to be "considered." The floating- 

point zero which makes up the last card shown in Figure 1. 2 indicates that there 
is no card-input measurement data for this case. 

Referring now to Figure 1. 3 under the heading "Parameter Estimation," we 
see eight columns of data. These are the measurements. Each row contains the date 
in the first two columns (YYMMDD , HHMMSS) , then three components each of position 
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and velocity (selenocentric, EE50). The data shown here were generated by the 
GTDS program and placed on unit 21 for processing by PEST. They represented 
osculating sets generated by perturbing certain elements of the JPL 15, 8 field 
truncated to 4, 4. It was the intent of this PEST run to determine which coef- 
ficients had been perturbed and by how much. The Method 7 initial conditions 
are printed, followed by the residuals in Method 7 elements. The residuals are 
in kilometers and radians, while time is in days from epoch. The estimated state 
increment tells how much change is to be made in the estimated parameters from 

^ —X ^ j, 

the residuals just printed. The penalty function printed out next is YQ Y for 
all of the observations. This is the function to be minimized by the (weighted) 
least squares procedure. The estimation process takes three trials to converge 
because the RMS state increment is greater than 10 for the first two trials. 

We may notice that the penalty function has reduced from . 5796 to . 7622 x 10 
after three iterations. The square-root information matrix is printed out in its 
triangular form (folded by format after 6 rows). This matrix provides infor- 
mation on the observability of the various parameters being estimated. The post- 
estimation gravitational coefficient values are printed out finally. These indicate 
that C(2, 0) , S(2, 1) , C(3, 0) , and S(3, 1) were perturbed (by comparing with the 
initially estimated field). 

Although the initial orbital element estimates are not printed, their incre- 
ments for the three trials add to zero (essentially), so we can assume that the 
residuals were entirely due to the four perturbed coefficients. The fact that the 
residuals (and penalty function) on the last trial were so near to zero indicates 
that no parameters other than those estimated -were perturbed. The measurement 
sensitivities for this process were computed only on the first trial, then "frozen" 
for the remaining trials without corrupting convergence appreciably. 

Figure 1. 4 shows input for a two-case run. No output is shown for this 
run. The first case is a simulation of measurements. Initial elements are taken 
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Figure 1. 4 Input for Simulation and Residual Print 


from locations 30-35 in INPUT common, then converted according to KMETH 
(1036) for propagation. The time step for measurements (location 460) must also 
be provided for this case. The first zero in column 6 causes escape from INPUTF. 
The zero in column 3 tells FILTER that a simulation is to be performed. The 
floating-point zero is to "end" measurement data input from cards and signal 
another case following. The next 16 zero is read by INPUTF at the beginning 
of case 2, causing a return to FILTER. Notice that more INPUT cards could 
have been read in before leaving INPUTF. FILTER reads the next card. The 
-1 in column 12 tells FILTER not to call WLS, but to compute and print residuals. 
The 2, 18, and 49 indicate that sensitivities of the measurements to C(2,0), 

C(2, 1) , and S(2, 1) are to be printed. The only necessary quantities on the 
measurement selector card for this case are M and IQ. The last zero, which 
signals "no more" measurement data cards, is omitted since this is the last case 
of the run. 
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n. PEST PROGRAM DESCRIPTION (Programmer's Guide) 

This section supplements Reference 2 by describing the PEST program 
insofar as it differs from MAESTRO. 

2. 1 Program Structure 

Figure 2. 1 is an overview of the structure of PEST and a flow chart of the 
main program or driver. The input and initialization functions for trajectory 
propagation are identical to those of MAESTRO. The primary difference is that 
PEST includes a path for MODE - 8 , the parameter estimation mode. A lunar 
lifetime mode (6) is retained from the MAESTRO capabilities in order to verify 
the characteristics of an estimated lunar field. For further elaboration of PEST's 
structure, see first the flow chart for FILTER (Section 2. 3) and then that for WLS. 

2. 2 Subroutine Cross-Reference and COMMON Blocks Descriptions 

Table 2. 1 lists the subroutines which make up the PEST program and defines 
their functions or purposes. It does not include those subroutines which are in com- 
mon with MAESTRO. 

Figure 2. 2 is a cross-reference map for the parameter estimation portion 
of PEST. For further details, see Figures 3. 1 and 3. 3 of Reference 2, 

Table 2.2 lists the COMMON blocks referred to by PEST subroutines. Those 
COMMON blocks marked by asterisks are found only in PEST — not in MAESTRO. 
The PEST -only COMMON blocks are described in Table 2. 3. The others are de- 
scribed in Reference 2. 
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PEST Main Program Flow Chart 


ENTER 



Figure 2. 1 PEST Program Overview 
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TABLE 2. 1 ' 


PEST SUBROUTINE DEFINITIONS 


MAIN 

ACCEL 

AVEQNS 

AVSTRT 

EQNS 

FIELD2 

FILTER 

GRAV 

MATMPY 

MODEL 

OSCAR 

OUTPUT 

PHOUSE 

PUTINV 

REDATA 

RKSEVN 

SHIMMY 

SOLVE 

SORDID 

SPNM 

SQRUT 

TIMEC 

TUMULT 

UPDATX 

WLS 


Initiates and drives PEST 
Determines spacecraft accelei’ation 

Averages derivatives of the elements over one revolution 

Averages osculating elements and variational equations 

Determines derivatives of the state for RKSEVN 

Evaluates central aspherical disturbing acceleration and 
its gradient 

Initiates the parameter estimation process 

Evaluates disturbing acceleration from noncentral bodies 

Matrix multiplication routine 

Controls state and sensitivity propagation 

Averages osculating elements 

Prints trajectory and sensitivity information 

Performs matrix reduction by Householder's algorithm 

Inverts a packed upper triangular matrix 

Initializes and reads measurement file 

Seventh-order Runge-Kutta integration routine 

Evaluates derivative for linear variational equations 

Solves weighted least squares equation for state increment 

Sorts and prints parameters to be estimated 

Evaluates Legendre polynomials for FIELD2 

Computes triangular square root of P-D symmetric matrix 

Controls logic for trajectory propagation 

Multiplies by a packed upper triangular matrix 

Updates state with the estimated increment 

Weighted least squares estimation routine 
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Figure 2.2 Parameter Estimation Cross-Reference Map 


REDATA — DATE, OSCAR , ORBIT, TRMN, MAPMPY 
MODEL — FOWARD , ORBIT, TRMN 


FILTER 


| REDATA 

MODEL 

SQRUT . 

PUTINV 

WLS — 

TUMULT 

PHOUSE 

SOLVE 

UPDATX 

OSCAR — M50LEQ, ORBIT, MVTRN, TRMN, PLANET, INTEG 
SORDID, ORBIT, TRMN, DATE, M50LEQ, MVTRN, CALEND 



TABLE 2.2 


COMMON BLOCKS 


COMMON 


PEST Subroutines Referring to this COMMON 


BLANK 

AVG 

CETBL2 

CETBL3 

*CNSDR 

CNTRL 

*COMIC 

CONST 

*DUMY 

*DXCOM 

FIELDM 

*FIELDX 

GRAVTY 

*HELMO 

*IMCOM 

INPUT 


INPUTS 

INTER 

INTVAR 

*INTVRX 


OUTPUT 
MAIN, AVEQNS 
MAIN 
MAIN 

FILTER, UPDATX 

MAIN, ACCEL, AVEQNS, AVSTRT, EQNS, FIELD2, GRAV, 
OUTPUT, RKSEVN, SHIMMY, TIMEC 

REDATA, UPDATX 

MAIN, AVEQNS, AVSTRT, EQNS, FIELD2, FILTER, GRAV, 
MODEL, OSCAR, OUTPUT, REDATA, SHIMMY, TIMEC, 
WLS 

MAIN, AVEQNS, EQNS, FIELD2, GRAV 

MODEL, PHOUSE, REDATA, SOLVE, UPDATX, WLS 

MAIN, FIELD2, FILTER, MODEL, SHIMMY, SORDH) 

UPDATX 

ACCEL, FIELD2, GRAV, OUTPUT, SHIMMY 

ACCEL, AVEQNS, EQNS, FIELD2, GRAV, OUTPUT 

FILTER, MODEL, REDATA 

FILTER, MODEL, REDATA, UPDATX, WLS 

MAIN, ACCEL, AVEQNS, AVSTRT, EQNS, FIELD2, FILTER, 
GRAV, MODEL, OSCAR, OUTPUT, REDATA, RKSEVN, 
SHIMMY, SORDID, TIMEC, UPDATX, WLS 

MAIN 

TIMEC 

MAIN, AVEQNS, AVSTRT, EQNS, FIELD2, GRAV, OSCAR, 
OUTPUT, RKSEVN, SHIMMY, TIMEC 

MAIN, AVEQNS, AVSTRT, EQNS, FILTER, MODEL, OSCAR, 
OUTPUT, REDATA, RKSEVN, SHIMMY, TIMEC, WLS 
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*IXCOM 

PERT 

PLNET 

SAVE 

STATE 

*YCOM 


MAIN, AVSTRT, FIELD2, FILTER, MODEL, OUTPUT, 
PHOUSE, REDATA, SHIMMY, SOLVE, SORDID, 

UPDATX, WLS 

MAIN, ACCEL, AVEQNS, EQNS, FIELD2, GRAV, SHIMMY, 
TIMEC 

MAIN, ACCEL, GRAV, OUTPUT 
MAIN, TIMEC 

MAIN, AVEQNS, AVSTRT, EQNS, FIELD2, FILTER, MODEL, 
OSCAR, OUTPUT, REDATA, SORDID, TIMEC, UPDATX 

FILTER, MODEL, REDATA, WLS 
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TABLE 2. 3 


Block 

COMIC 

DUMY 

DXCOM 

FIELDX 

HELMO 

IMCOM 

INTVRX 

IXCOM 


PEST COMMON BLOCK DESCRIPTIONS 


Length 

(DP words) Use and Make-up 


6 

27 


5850 


309 


12 ' 


4 


1212.5 


54 


Used to transfer estimated initial elements 
ELMI(6) Estimated initial elements 

To communicate elements, transformations, rates 
DUM(6) Dummy element vector 

BUMROT(3,3) Equatorial to orbital transformation 
ELMD<12) Storage for elements, sines, and 

cosines 

Communicates crucial estimation information 
DX(100) Estimated state increment ^ ^ 

HQY(IOO) Weighted state residual, H^Q S Y 
QH(5650) Square-root information matrix and 
weighted sensitivities 

transfers variational equation derivatives 
DVDX2(3, 3) Gradient of the disturbing acceleration 
DFDCX(3, 100) Explicit sensitivity of gravitational 
disturbing acceleration to estimated 
parameters 

Transfers estimated elements and secant increments 
ELMO(6) Current estimated orbital elements 

EODINK(6) Fraction of initial elements used in 
computing secant partials 

Communicates measurement indicators, gain, tolerance 
IM(6) Measurement indicators 

IGAIN Reciprocal gain for including state increment 
LTOL Iteration convergence tolerance (negative -tog-j^) 

Communicates measurement sensitivities and rates 
H(60G) Measurement sensitivity, by/bx 

QATES(606) Rates of H for variational equations 
NQ Number of variational equations 

Communicates estimation control indicators 
IX(100) Identification array for estimated parameters 

N Number of parameters to be estimated 

Ml Number of elements observed at one time 

M Number of observations (sets of size Ml) 

NQH Number of elements in the square-root 

information matrix N*(N+l)/2 


26 



YHAT 


IQ Weighting treatment key 

ITMAX Maximum iterations allowed 

ID Observation counter (1 ^ ID SM) 

IT Iteration counter (1 ^ IT ^ ITMAX) 

50 Communicates time, measurement, estimate, and weight 

TI Current time 

YHAT(6) Estimated elements 

TF Next measurement time 

Y(6) Measured elements 

Q(36) Measurement weighting matrix 
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2. 3 Subroutine Descriptions 


Documentation is included here for those subroutines which are uniquely 
required by PEST. PEST requires many of the subroutines of MAESTRO, but 
since these are described in detail in Reference 2, only differences between 
MAESTRO and PEST versions of the common routines will be discussed. Table 
2.4 categorizes the subroutines. 


TABLE 2.4 

PEST/MAESTRO SUBROUTINES 


PEST Subroutines Modified MAESTRO Subroutines 


FILTER 

MODEL 

OSCAR 

PHOUSE 

PUTINV 

REDATA 

SHIMMY 

SOLVE 

SORDID 

SQRUT 

TUMULT 

UPDATX 

WLS 


MAIN 

ACCEL 

AVEQNS 

AVSTRT 

EQNS 

FIELD2 

GRAV 

MATMPY 

OUTPUT 

RKSEVN 

SPNM 

TIMEC 


The primary modification to the MAESTRO subroutines is to increase the 
size of INTVRX common from 13 to 1213 to permit calculation of 100 sets of 
linear variational equations. Other differences are noted below. 

MAIN FILTER is called when MODE = 8 

ACCEL Contains FIELDX common and initializes DVDX2 array for 

linear variational equations (LVE) computation 
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AVEQNS 

AVSTRT 

EQNS 

FIELD2 

GRAV 

MATMPY 

OUTPUT 

RKSEVN 

SPNM 

TIMEC 


Averages linear variational equations as well as elements 

Averages initial LVE as well as initial elements 

Calls SHIMMY in computation of LVE 

Computes acceleration sensitivities to gravitational co- 
efficients for LVE 

Computes gradient of non-central-body acceleration for LVE 
Includes a C = AB T call (MATMPT) 

Includes extra output for LVE 

Has larger internal storage for INTVRX variables 
Generates Legendre polynomials a row-at-a-time for FIELD2 
Only differences are for larger INTVTtX common 



SUBROUTINE FILTER 


Calling Sequence: 


CALL FILTER 


Purpose: FILTER initializes the estimation process and 

controls (1) simulated data generation and (2) 
nonestimation sensitivity and residual calculations 

Common Blocks Required: CNSDR, CONST, FIELDM, HELMO, IMCOM, 

INPUT, INTVRX, IXCOM, STATE, YCOM 


Subroutines Called: CALEND, DATE, MODEL, MVTRN, M50LEQ, 

ORBIT, OSCAR, REDATA, SORDID, TRMN, 
WLS 


Input/ Output 



I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

I 

PI, PI2 

1 

CONST(2, 3) 

77, 277 

I 

GM 

12 

CONST<5) 

Planetar}' gravitational constants 

I 

ZONL 

288 

FIELDM(l) 

Lunar field coefficients 

0 

ELMO 

6 

HELMO(l) 

Estimated orbital elements 

o 

H 

606 

INTVRX(l) 

Measurement sensitivities 

m 

T FINAL 

1 

INPUT (4) 

Final time 

u 

cov 

21 

INPUT(56) 

Input weighting matrix 

I 

A(197) 

2 

j INPUT(197) 

Solar pressure coefficients 

i 

STEP 

1 

1 

t 

| INPUT(460) 

i 

Simulated measurement time 
increment 

0 


1 

f 

INPUT(1032) 

Integration stop type key 

I 

KMETH 

1 

INPUT(1036) 

Integration method key 
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I/O 


Symbolic 

Name 


Program 

Dimension 


Common 

Block 


Dimension 


a 

IX 

100 

IXCOM(l) 

Estimated state identifiers 

0 

N " 

1 

IXCOM(lOl) 

Number of state components 

0 

MI 

1 

IXCOM(102) 

Number of observables at each 
measurement 

0 

Y 

6 

YCOM{9) 

Propagated elements 

0 

Q 

36 

YCOM(15) 

Measurement weighting matrix 

o 

ICNSDR 

20 j 

CNSDR(l) 

"Considered" parameter identifiers 

■ 

DJO 

1 

INPUT(46) 

Initial Julian date 

I 

KVAR 

1 1 

INPUT(1096) 

! 

Key for sensitivity generation 
method 

I 

KAVST 

— 

1 | INPUT(1098) 

Key for averaging elements, 
setting I. C. 's 

0 

NQ 

1 

INTVRX(1213) 

Number of sensitivity elements 

I 

IM 

6 

IMCOM(l) 

Measurement selector array 

I 

IGAIN 

1 

IMCOM(7) 

• 

State increment inclusion number 

• 

B 

LTOL 

1 

IMCOM(8) 

Estimation iteration tolerance 
number 

o 

ID 

1 

IXCOM(107) 

Measurement counter 
(set for REDATA) 

o 

IT 

1 

IXCOM(IOS) 

Iteration counter 

0 

EJO 

1 

STATE (2 6) 

Ephemeris time epoch 

0 

M 

1 

IXCOM(103) 

Number of measurements to be 
processed 

o 

NQH 

1 

IXCOM(104) 

Number of elements in SRIM 

o 

IQ 

1 

IXCOM(105) 

Weighting type key 

0 

IT MAX 

1 

IXCOM(106) 

Maximum iterations in WLS 

I 

X 

6 

STATE (1) 

Initial estimated Cartesian state 

0 

TI 

1 

YCOM(l) 

: Previous measurement time 







O 


TF 


YCOM(8) 


Current measurement time 















































































Discussion 

FILTER performs three basic functions 

1. It initializes for WLS, the estimator 

2. It contains logic for generating a simulated data file. 

3. It contains logic for reading the data file, printing measurements, 
residuals and sensitivities. 

The IX array is read in from the input unit. It is this array which determines which 
function FILTER will perform: 

IX(1) < 0 , N= 0 , generate simulated data file 
IX{1)>0, IX(N + 1) 0, call the estimator, WLS 

IX(1) > 0 , EX(N + 1) < 0, read data file, compute and print 
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FILTER FLOW CHART 
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Return 
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* 


Return 
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SUBROUTINE MODEL 


Calling Sequence: 


CALL MODEL 


Purpose: MODEL computes the estimated measurement 

(orbital elements) and secant partials of the 
measurement with respect to state at a specified 
time. 

Common Blocks Required: CONST, DXCOM, FIELDM, HELMO, IMCOM, 

INPUT, INTVRX, IXCOM, STATE, YCOM 

Subroutines Required: FOWARD, ORBIT, TRMN 

Input/ Output 


I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

I 

GM 

12 

CONST{5) 

Gravitational constant array (GM) 

o 

QH 

5650 

DXCOM(201) 

: 

Measurement sensitivities to state 

I/O 

ZONL 

1 

16,18 

FIELDM(l) 

Lunar field harmonic 
coefficients 

I/O j ELMO 

6 

HELMO(l) 

Estimated elements 

I 

EODINK 

6 

HELMO(7) 

Initial element increments for 
partials 

I/O 

H 



606 

INTVRX(l) 

Measurement sensitivities to state 



i— 1 

INTVRX(1213) 

Sensitivity computation key 

0 

T FINAL 

1 



INPUT(4) 

Integration stop time 

1/0 

A123(198) 

1 

INPUT(198) 

I Solar radiation pressure 
coefficients 

I 

KMETH 

1 

INPUT(1036) 

Integration method type 

I 

KVAR 

1 

INPUT(1096) 

Sensitivity computation type » 

key 

m 

IX 

N 

IXCOM(l) 

State component identifier 
array 

i 

N 

1 

IXCOM(lOl) 

1 

Number of state components 
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I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

I 

NQH 

1 

IXCOM(104) 

Number of square-root matrix 
elements 

0 

X 

6 

STATE(l) 

Cartesian state of spacecraft 

I 

X<36) 

1 

STATE(36) 

Solar pressure coefficient 

0 

T 

1 

STATE(IO) 

Time 

0 

ELM 

6 

STATE(14) 

Orbital elements of 
spacecraft 

I 

TI 

1 

YCOM(l) 

Entry time 

0 

YHAT 

6 

YCOM(2) 

Estimated measurement 
(elements) 

■ 

TF 

1 

YCOM(8) 

Exit or measurement time 

I 

IM 

6 

IMCOM(l) 

Measurement identifiers 


Discussion. 

MODEL calls FOWARD to propagate the estimated orbital elements, ELMO, from TI 
to TF. Secant partials of the estimated measurement, YHAT = ELMO, are computed 
in H and in QH beginning with QH(NQH+1) . The reason for double storage is that the 
QH array is destroyed by the estimation process in subroutine PHOUSE, while H is 
used to store initial conditions for integration between measurements. That is 


H.. = S 

i] 


bYHAT.(T) 

5CL 

J 


for gravitational or solar pressure coefficients, C. 


bYHAT.(T) 

■ ■_ T , '• — for initial condition states, ELMO.(O) 
OELMO^(O) j v ' 


r 


and is computed from 


YHAT(TF, S + A) - YHAT(TF, S) YHAT(TF, S + A) - ELMO(TF) 
' } A A 
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where S Ls the estimated state and A is a small change in S . Rather than 

A A 

carrying along YHAT(T, S+A) from measurement to measurement, YHAT(TI, S + A) 
is reconstructed from H to initiate integration over the next interval. 

YHAT(TI,S + A) = ELMO(TI) + H(TI) * A 

If KVAR is nonzero, linear variational equations replace secant partials as 
the means for computing measurement sensitivities. These sensitivities reside in 
H (INTVRX) upon return from FOWARD. 

The QH array (DXCOM) is upper-loaded from H according to the measure- 
ment identification array, IM. H is viewed as a 6xN array, while QH is viewed 
as an Mlx6 array, although both are programmed as single-dimensioned arrays. 
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MODEL Flow Chart 



Load H-QEM 


< 


RETURN 


3 












SUBROUTINE OSCAR 


Calling Sequence: CALL OSCAR(KIN, TI, El, TO, EO) 

Purpose: OSCAR converts osculating input (elements 

or Cartesian state) to averaged elements. 


Common Blocks Required: CONST, INPUT, INTVAR, INTVRX, STATE 

Subroutines Called: INTEG, MVTRN, M50LEQ, ORBIT, PLANET 

Input/ Output 


I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

I 

KIN 

1 

Call LIST 

- 

Input type key 

(8: Cartesian; other: elements) 

I 

TI 

1 

Call LIST 

Osculating reference time 
from epoch 

I 

El 

6 

Call LIST 

Input osculating elements 

O 

TO 

1 

Call LIST 

Output reference time from 
epoch 

0 

EO 

6 

Call LIST 

Averaged elements output 

I 

RAD, PI2 

1 

CONST(l) ! 

Angular constants 

I 

GM 

12 

CONST{5) 

Gravitational constants 

I 

DJO 

1 

INPUT(46) 

Julian date at epoch 

I 

KMETH 

1 

INPUT(1036) 

Trajectory propagation method 
key 

. 

1 

I 

] 

XX 

6 

STATE(l) 

Input osculating Cartesian 
state 

I/O 

X, H 

1 

INTVAR(l) 

Integration time and step 

I/O 

Z 

6 

INTVRX(l) 

Intermediate osculating elements 
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Discussion 


OSCAR is very similar to AVSTRT, but may be called at times other than epoch 
to average osculating sets other than the estimate. It is used to average "measured ' 1 
osculating sets for use with averaged propagation methods in estimation. The os- 
culating elements are integrated forward over one mean orbital period (found by 
iteration). The integral is then divided by that period to render the averaged ele- 
ments, which are taken to be referred to the center of the integration interval. 



OSCAR 



Initialize for averaging. 
Save integration variables. 


Changed initial conditions, 
so call INTEG(-l) 



100 ] — * 


Call INTEG to take one R-K step. 
Compute period of average orbit, 
compare it with the current time 


7 ^ 


Time is greater than t+TOL. 
Compute step size to make T = r 
using Newton-Raphson formula. 
Re-initialize at previous time, 
continue 


No 


No 

Save time and 
_elem£nts 


T> T- 

TOL? 

> 

Yes 

r 


|T-r j<TOL? 


Yes 



Divide elements by T to average. 
Set output time to T/2. 

Restore saved parameters. 





RETURN 


3 


i 
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SUBROUTINE PHOUSE 


Calling Sequence: 

CALL PHOUSE 

Purpose: 

PHOUSE reduces a partly-upper-triangular 
matrix to purely-upper-triangular, using 
Householder's Reduction Algorithm. 

Common Blocks Required: 

DXCOM, IXCOM 

Subroutines Called: 

None 

Input / Output 



r~ . . . 1 

I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

I 

N 

1 

IXCOM(lOl) 

Columns of input and output 
matrix 

I 

MI 

1 

IXCOM<102) 

Rows in excess of N 

I 

NQH 

1 

IXCOM(104) 

N * (N+l)/2, non-zero 
elements of N x N 

— 

I/O 

QH, WT 

NQH+N*M.l 

DXCOM(l) 

Square-root information matrix 
(packed upper-triangular) 


Limitations: QH is singly-dimensioned as follows 


QH " (W 11- W 12' *22’ W 13' W 23' ••• W m ' h ll- h 21’ ••• h mX’ h 12’ h 22’ 


h , 0 — 
m,n 

(NQH+N*Ml)-th 

element 


I 

NQH-th 

element 


The input QH array is replaced by the output array. Ml must not exceed 6 and NQH+ 
N*M1 must not exceed 105 unless program dimensions are changed. 

Method 

PHOUSE is a special-purpose implementation of Householder's Reduction Algorithm (1,2) 
for reducing a matrix to upper-triangular form. It is used in the weighted least squares 
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T -1 

process to accumulate the information matrix, HQ H, in upper-triangular square- 
root form. In this discussion, we can let Q 1 = I without loss of generality. 

Let us assume that at some point in a WLS process we have an upper-triangular SRIM 
(square-root information matrix), W, which is N.x N. That is, W contains N x (N+l)/2 
potentially non-zero elements above the diagonal and zeros below. W also contains all 
of the information from k-1 measurements. 


T 

W W 
k-1 k-1 


k-1 

Z 

i = 1 


T 

H. H. 

l i 


( 1 ) 


In the equation, H. is the sensitivity of measurement to state at the i-th measurement. 
We now consider what happens to W when we add another measurement, whose 
sensitivity matrix is H^. The new SRIM, W^, must obey 


W. T W, = W t T W, , 
k k k-1 k-1 


+ 


( 2 ) 


which could be written 


w k w k - 

T Tl 

W H 

k-1 k 

* 

i 

1 


J 

|> J 

(nxn)(nxn) 

nx(n+m ) 

(n+m)xn 


We may look at the reduction problem then as the problem of finding an orthogonal 
matrix, R, such that 


"w. 

R 

W. , ’ 

k 


k-1 

0 


H k 


because 


* 

T 




T 


w. 


w, 


w, . 

_T 

w, , 

k 


k 

— 

k-1 

R R 

k-1 

0 


0 


> _ 


_«k 


(4) 


(5) 


is identical to equation (3) for orthogonal R. Let us now simplify the notation by 
dropping the subscripts and renaming so that equation (4) becomes (6). 


W = RA 


( 6 ) 
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Householder's choice of R is the product of a sequence of reflection matrices (which 
are symmetric and orthogonal). 

(Introduction of new indices is unavoidable) 

R^) defined by 


R 


« - I - /3,w‘ k >W< k > T 


( 8 ) 


where I is the (n+m) identity, W ^ is a column vector to be developed, and 


*k = 2/ |w (k) | • 


(9) 


The derivation of and R^ \s inductive. If A of equation (6) is partitioned by 
columns, a., as in equation (10) 

A = [ a l a 2”-' a n] (10) 

we note from equation (2) that the first element of W must be - a^. We define 

“l = | a i I Sign 


If we choose R^ to zero the subdiagonal elements of where is defined 


by equation (12), 

w< k > - H«R<*- 1 >....H< 1 >A B rI "'/' 11 
we may solve for (which turns out to be the k-th column of W^). 


(12) 


>(1). 


ri^M) 


R Vi/ a 1 =' a x (3 ± W w W v ' ^ = a ± e t 


(13) 


In equation (13), e is the first unit basis vector. The unit basis vector e. has 1 
1 \ 1 


(1) 

as its i-th component and zeroes elsewhere. We may solve (13) for W v , 


W 


-< 1 ) 


a i •*’ “i °i 

w ■ 


(14) 


(k) 

but we do not know how to compute the denominator of (14). The length of W v ’ in 

equation (8) is arbitrary, however, due to the definition of (3 , so we can simply 

T (1) 

define its length to be such that a^ W' ' = 1. Then we have equation (15). 


45 



(15) 


w<1) = a i + a i e i ■ 

Having now solved for W^, we could form from equation (12) by making the 
association, « A. The first column of W^ is zero, except for the first element, 


which is equal to - a , We could now repeat the process, equation (12), by looking for 
which would zero the subdiagonal elements of the second column of without 
disturbing the first column. We can preserve the first column of W^ under the 

' to 

( 1 ) 


(2) (2) 

transformation, R v , by setting the first element of W w to zero. Then we have 


R 


^ W^ = W^ - )3 

2 



= W' 


The construction of is otherwise identical to that of W^. 


we can arrive at the following set of equations. 

■ 

sign 


ot =\ 

k 


' m 


kk 


(16) 


Proceeding inductively, 


(17) 


0k 


ot (a. + ) 

k^ k kk' 


(18) 


w: 


(k) 


i <k 


,( k ) . 


a + w v 

k kk ’ 


W 


(k) 

ik 


i = 


i = 


k 

k+1, n+m 


(19) 


(k) 




0 

1 

( k ) T w! k) 


i < k 
i = k 

i = k+1, n+m 


W. (k) in equations (20) and(21) is the i-th column of W^. 


^■(k +1 ) _ ^y(k) _ Y^k) ^(k) 

i i i 


i = k+1, n+m 


( 20 ) 


( 21 ) 


Congratulations to you if you have followed the development to this point. It was not 
easy to write or type it, either. 
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In coding this routine, some advantage has been taken of the fact that the initial 

matrix is partly upper triangular at entry. That is, W 1 of equation (4) is upper- 

(k) k_1 

triangular to begin with, so the dimension of W^' * in equation {19) needs only 
to be m+1 and the inner product of equation (20) needs to sum only m multiplications. 
These advantages can be significant if n is large. 
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SUBROUTINE PUTINV 


Calling Sequence: 

Purpose: 

Common Blocks Required: 
Subroutines Called: 


CALL PUTINV (B,C,N) 

PUTINV inverts a packed upper triangular 
matrix. 

None 

None 


Input/Output 


I/O 

Symbolic 

Name 

Program 

Dimension 

Communication 

Definition 

I 

B 

N(N+l)/2 

h“ — 

Call List 

Input matrix (P.U.T,) 

0 

C 

N(N+l)/2 

Call List 

Output matrix (inverse of 
B, P.U.T.) 

I 

N 

I 

Call List 

Rows or columns of upper 
triangular matrix 


Note: B and C are single -dimensioned arrays containing the potentially non-zero 
elements of an upper triangular matrix. 


B (b u , b 12 , b 22 , b 13> 


, b ) 
nn 
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Method 


The input matrix, B, and its inverse,. C, are treated as upper triangular arrays. 


same for C. 



b u 

b 12 

b 13 

etc. 

B = 

0 

b 22 

b 23 

etc. 


0 

o 

etc. 

b 33 

etc. 


They must obey CB = I, which is equivalent to the following set of equations. 


CB =1 


Solution 


°n b n = 1 

C 11 =1//b ll 

C 22 b 22 “ 1 

C 22 =1 ^ b 22 

C 11 b 12 + °12 b 22 “ ° 

C 12 = ~ C 22 °1 

°33 b 33 = 1 

C 33 ~ 1//b 33 

C 11 b 13 + C 12 b 23 + 'l3 b 33 ° ° 

C 13 ” _C 33^ C 1 

c 22 b 23 + c 23 b 33 ~ ° 

°23 = ~ C 33 (C : 

etc. 

etc. 

n 

c kk = 1/b kk 

y' c..b =6 

jk ik 

1 

c jk c kk( a 


"12 23' 


k-1 


*=i 


B and C are treated as single - dimensioned arrays, so the programming requires 
calculation of the appropriate single indices which correspond to the double indices 
used in the development. 
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SUBROUTINE REDATA 


CALL REDATA 

REDATA reads the data file and loads 
measurements into common. It -also 
initializes certain variables at the be- 
ginning of the case and at the beginning 
of each iteration. 

Common Blocks Required: COMIC, CONST, DXCOM, HELMO, IMCOM, 

INPUT, INTVRX, IXCOM, STATE, YCOM 

Subroutines Called: DATE, MATMPT, MATMPY, ORBIT, 

OSCAR, TRMN 

Input/Output 


Calling Sequence: 
Purpose: 



I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

I/O 

EL MI 

6 

C0MIC{1) 

Estimated initial elements 

0 

QH 

5650 

DXCOM(201) 

Initial sensitivity matrix for 
WLS 

I/O 

ELMO 

6 

HELMO(l) 

Initial elements for integrator 

o 

EODINK 

6 

TT _ T ■ Initial element increments for 

HELM0 < 7 > I secant partials 

o 

H 

606 

T\m 7 i>v/i\ i Initial sensitivity matrix for 

INTVRX(l) { M0DEL 

m 

N 

1 

t 

IXCOM{101) [ Number of states 

I 

Ml 

1 

< Number of observables at each 

iAL^UiVl(lUZ) f 

f measurement 

o 

| M 

! 1 ! 

IXCOM(104) 

Number of measurements 

I 

NQH 

1 

IXCOM(105) 

Number of elements in SRIM, 
N{N+l)/2 

I 

IQ I 1 

IXCOM(106) 

Weighting flags 

I 

ID 

1 

IXCOM(IOS) 

Measurement counter 

I 

IT 

1 

IXCOM{109) 

Iteration counter 

0 

TI 

1 

YCOM(l) 

Most recent measurement 
time 
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I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

o 

YHAT 

- 6 

YC0M(2) 

Estimated measurement vector 

I 

GM 

12 

CONST(5) 

Gravitational constants 

I 

DJO 

1 

INPUT(46) 

Julian date at epoch (UMT) 

I 

KMETH 

1 

INPUT(1036) 

Trajectory propagation method 
key 

I 

KAVST 

1 

INPUT (1098) 

! 

Averaging indicator, 
start-up key 

0 

EJO 

1 

1 STATE(36) 

1 

Ephemeris date at epoch 

I/O 

X 

6 


Cartesian state 

I/O 

TF 

1 

j YCOM(8) 

Measurement time 

0 

Y 

6 

YCOM{9) 

Measurement vector 

o 

Q 

36 

YCOM(15) 

Measurement weighting matrix 

I 

IM 

6 

IMCOM(l) 

Measurement identifiers 


Discussion 

Subroutine REDATA perforins two functions: (1) loading the measurement from the 
data file and (2) initialization of certain program variables. The data file (logical 22) 
is read once on each entry with an unformatted binary read, expecting a time and six 
elements on each record. When an end-of-file is encountered, M , the number to signal 
"no more data — estimate the state now 1 ' is sent to WLS. The initialization functions are: . 

1. Load the estimated initial element vector with elements developed 
from inputs at the beginning of the case. 

2. Load the integrator's initial element vector, ELMO, and the estimated 
initial measurement vector, YHAT, with the estimated initial elements 
at the beginning of each iteration or pass through the data. 

3. Rewind the logical unit on which the data reside at the beginning of 
each iteration. 
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4. At the beginning of each iteration, set the measurement sensitivity 
matrix (H and part of QH) to zero — or to one in the case of the 
sensitivity of an initial element measurement to an initial element 
state. The secant increment is also set when an initial element 
is to be estimated. The increment is 1/100,000 of the element 1 s 
initial value. 

Initialization Logic (see flow chart) 

Input measurement data is assumed to reside on unit 21 in the form of calendar date 
plus Cartesian state (EE50). The state may be regarded as osculating or averaged. 

It will be converted to processing form depending upon the propagation method, KMETH, 
and the averaging key, KAVST. The time from epoch and the converted state (elements 
of the KMETH type) are written onto logical unit 22. If KMETH is 4 or 7 , the 
measurements need be converted only once, but if KMETH is 5 or 8 (averaging), 
the measurements will be re-averaged at each iteration with the current estimated 
lunar gravitational field. 

The weighting matrix, Q, used in weighted least squares estimation under the IQ = 3 
option can be converted from Cartesian EE50 variables to method KMETH variables 
by REDATA. It is expected that COVIN will be read in from logical 21 if IQ = 3, and 
that COVIN will be the upper (or lower) half of a covariance matrix in the Cartesian 
state. This is loaded into a filled-out, symmetrical Q and converted to elements by 



The required partial derivatives are computed by the secant method. Q is written 

E 

onto logical 22 for potential use by subroutines WLS and SQRUT. 

The number of measurements written onto 22 will be the smaller of M 
(input maximum number) and the number found on logical 21. 
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REDATA Flowchart (I) 



Initialize estimated 
initial elements 
from input 

Load initial elements 
from (updated) estimate 


Load measurement 
from data file 


Set flag to indicate 
end of data 
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REDATA FLOW CHART (H) 
(Initialization Logic) 
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SUBROUTINE SHIMMY 


Calling Sequence: 
Purpose: 


Common Blocks Required: 


CALL SHIMMY (EM, D, ELMR, R) 

SHIMMY calculates the rates in the linear 
variational equations for the variables p, 
ecos (u$, e sin (a), to + f, i, Cl 

CONST, CNTRL, FIELDM, FIELDX, INPUT, INTVAR, 
INTVRX, IXCOM, PERT. 


Subroutines required: MATMPY, ROTATE. 


Input / Output 


M j 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

i 

D 

3,3 

Calling 

Operand 

Transformation matrix from 
M50 to position vector 
coordinates. 

i 

DFDCX 

3,100 

FIELDX (10) 

Partials of F w.r.t. 
harmonic coefficients 

i 

DVDX2 

3,3 

FIELDX (1) 

Gradient matrix 

i 

ELMR 

12 

Calling 

Operand 

Orbital elements with sines 
and cosines 

i 

EM 

6,3 

C ailing 
Operand 

Matrix which relates the 
disturbing acceleration to 
the orbital element rates 

i 

GM 

12 

CONST (5) 

Gravitational constants 

i 

NB 

1 

IXCOM (101) 

Number of coefficients to be 
studied 

0 

QATES 

600 

i 

INTVRX (601) 

The rates of the linear 
variational equations 

a 

R 

i 

Calling Operand 

Distance to S/C 

■ 

RCART 

! 

3 

PERT (1) 

Disturbing acceleration of 
S/C in mean of 50 
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Input / Output (continued) 


1 

RSW 

3 

PERT (4) 

Disturbing acceleration of 
S/C in orbital 
coordinates 

I 

SELNEQ 

3,3 

FIELDM (289) i 

j 

Transformation matrix from 
SELNEQ=Q T to inertial 
selenographic coordinates. 

I 

Z 

600 

INTVRX (l) 

1 

Linear variational 
equation variables 


Theory: The equations for the time rates of the orbital elements are 

represented by: 

E - F <E,C, r (E), t), 

where r (E) represents an implicit dependence upon the position/ velocity vector. 
These equations may be differentiated to yield equations which represent the 
effects of small variations in the gravity coefficients upon the subsequent history 
of the orbital elements} 

ie. jL /A _d_ dE_ J/3>F\ JbA 

be \E J dt 6C j \dE ) \bfj j dC \bCj , 

where the normal brackets ( ) represent explicit differentiation. The equations 
of motion can be written in matrix form as follows: 


t 


' 


1 


R 


- — 





E 

6x1 


N 
6x 1 


M(E) 

6x3 


S 

05 

= [N] +[M(E)] 

D(E) 

3x3 


Q(t) 

3x3 


Tv" 

3x1 

* * 





3x1 


— 


- — 




where D(E) is the transformation matrix from the inertial frame to the radial, cir- 
cumferential, orbit normal orbital coordinate system, Q(t) is the rotation from seleno 
graphic to inertial coordinates, and 7 V is the gradient of the selenocentric potential 
function. The vector N contains the unperturbed (Keplerian) part of the rate of the 
true anomaly. The 6x6 matrix in curly brackets is evaluated in two parts made up 
of the explicit partial derivative matrix. 
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f ^ 

d ( P»h, k, u, i, 6) 
p, h, k, u, i, ) 


and the product 


(») (») ■ 


. where 


(I?)- HHMNK] 

If “ [ dT (s)] - (|e [ dT ](I) + [ DT Hi 

3 F ' 

The last three columns of rs, contain zeros because it is assumed that no 

dr ^ 

velocity-dependent forces are acting. Therefore, the last three rows of -^7- 

are not required and set to zero. The 3X3 matrix J is the Jacobian matrix 

£> 


of the disturbing potential with respect to the body-fixed cartesian coordinates. 
Finally, the matrix gives the explicit partials of the planetary equations 

with respect to the gravity coefficients. 



Description: The equations discussed above are of the form Z = AZ + B. 

SHIMMY calculates the matrices A and B, and then uses the input value of Z to 
calculate the new value of (QATES). The matrices Z and QATES are dimen-, 
sioned 6 X 100, and represent the quantities dE/dC and dE/dC respectively, 
where each column of the matrices corresponds to a different harmonic coefficient 
of the field. The input quantity NB determines how many coefficients are to be 
studied, with a maximum of one hundred. 


The matrix B represents the quantity 
written as 

(* 

where M 


lty (If)' dl 


discussed above, which may be 


UKd! 

r Q 1 T 52V 1 

) 1 J _ . 

J a c 

f — — — u» 

a 6 X 3 matrix which represen 


s the relationship between the per- 


turbing acceleration and the rates of the orbital elements. The matrix is calculated 


[d] & Q 


are 


in EQNS, and stored in EM. The transformation matrices 
discussed in the theory section and are stored in the 3x3 matrices D and SELNEQ 
respectively. Both of these arrays are calculated outside of SHIMMY. The matrix 


d VV 

a C 


is calculated in FIELD2 and stored in DFDCX. For further discussion 
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of this matrix and the order in which the elements are stored, please see 
subroutine FIELD2. 


The matrix A represents the quantity $F/ 3 E, and may be expressed as 


>(tt)S - 


Both of these quantities represent the product of 


several matrices, as indicated in the theory section above. In addition to the 
matrices already used for the calculation of B, the quantities discussed below 
must be available for the calculation of A. The matrix D2VDX2, calculated in 
FIELD2, contains the gradient of the force field and is used in the calculation of 
feF/gr). The vector RSW represents the disturbing acceleration on the S/C 
in orbital coordinates, and the vector RCART represents the same disturbing 
acceleration in inertial coordinates. The matrix , ie ( , * * U ’ 2 ) , 

\c E / \SP. h > k > i.n / 

is calculated column by column internally through successive matrix multiplica- 


ix (S) • “ 


lions and stored in DFDEX. The matrix and vector multiplications are done in 
subroutine MATMPY and ROTATE. Once these calculations are completed, 
the vector Z is inserted into the equation and the vector QATES is output. 
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SHIMMY FLOW CHART 
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SUBROUTINE SOLVE 


Calling Sequence: CALL SOLVE 

T 

Purpose: SOLVE solves the equation W WX = Y for 

X, given the vector, Y, and the upper-tri- 
angular matrix, W. 

Common Blocks Required: DXCOM, IXCOM 

Subroutines Called: None 


Input / Output 


r - . . i 

I/O 

Symbolic 

Name 

. Program 
Dimension 

Common 

Block 

Definition 

i 

QH 

NQH, max 105 

DXCOM(l) 

Square-root information matrix, W 

i 

KQY 

N, max 10 

DXCOM(106) 

Weighted residual vector, Y 

o 

DX 

N, max 10 

DXCOM(llG) 

Incremental state vector, X 

I 

N 

1 

IXCOM(lOl) 

i 

Number of states 


Note: The square-root information matrix, QH, is the singly-dimensioned representa- 
tion of a doubly -dimensioned packed upper-triangular matrix, W. 


QH = ( W u , W 12 , W 22 , W 13 , W 23> ..., W nn , + unused words) 

Method 

The weighted least squares estimation equation can be stated in the form 
(W T W) X = Y 

where W is the upper-triangular square-root information matrix, X is the estimated 

state increment to be solved for and Y is the weighted residual vector. An upper 

triangular matrix is a square matrix which has only zeroes below its main diagonal, 

T 

One obvious solution to the stated problem would be to form W W, invert it, and 
multiply Y by it to Obtain X. This subroutine uses another method, capitalizing on 
the fact that W is upper-triangular. We first solve the problem 
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T 

W B = Y 

for the dummy vector, B, by forward substitution. We then solve 
WX = B 

by backward substitution to obtain X. 


Forward substitution: 
T 

W B = Y 


Solution 


w ii b i = y i 
W 12 b l + "22 b 2 = y 2 


b i = y i /w u 

b 2 ' <y 2 - W 12 b l )/W 22 


) w..b. = y 
i=l 1J 1 


Backward Substitution: 
WX = B 


w x = b 
nn n n 


w , x + w x = b 
n-l,n-l n-1 n-l,n n n-1 



Solution 


x = b /w 
n n nn 


x = (b - w x )/w 
n-1 n-1 n-l,n n n-1, n-1 


n 


Iw..X. 

. . 1 

1 =] 



X. 



w,.x.)/w.. 

1 JJ 


Since the zeroes of W do not enter the solution, the potentially non-zero elements 
of W are stored in a singly-dimensioned array. The coding synthesizes the 
doubly-dimensioned indices to effect the above procedure as shown. 
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SUBROUTINE SORDID 


Calling Sequence: CALL SORDID 

Purpose: SORDID sorts the parameter identification 

array, IX, for FIELD2's required order 
and prints it out. 

Common Blocks Required: FIELDM, INPUT, IXCOM, STATE 

Subroutines Called: None 

Input/ Output 


I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

I 

ZONL 

16 

FIELDM(l) 

Zonal harmonic coefficients 

. I 

TSRL 

16, 17 

FIELDM(17) 

Longitudinal harmonic 
coefficients 

I 

NMOD 

1 

FIELDM(298) 

Maximum modeled N 

I 

MMOD 

1 

FIELDM(299) 

Maximum modeled M 

H- 1 

KMETH 

1 

INPUT(1036) 

Trajectory propagation method 
key 

I/O 

IX 

100 

IXCOM(l) 

Parameter identification array 

I 

N 

1 

IXCOM(lOl) 

Nonzero elements of IX 

I 

SPC 

1 

INPUT(198) 

Solar pressure 
i coefficient (280) 

I 

SPC 

1 

STATFnm 1 Solar P ressure 

STATE(36) | coeffiolent (278) 
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Discussion 


The IX array specifies the parameters to be estimated according to the code: 

1 - 272 Gravitational parameters 
IX(I) = < 273 -* 278 Initial elements 

[ 279, 280 Solar pressure coefficients 

The code for the gravitational parameters, C and S , is that IX(I) shall be 

ran nm 

the element number of the ZONL, TSRL array shown in the following figure. SORDID 
checks that the input values for IX are valid. That is, IX numbers must lie between 
1 and 280 and must not indicate coefficients outside the prescribed model (NMOD, 
MMOD). Furthermore, subroutine FIELD2 expects a particular order for the selected 
coefficients in order to efficiently evaluate sensitivities when the option for linear 
variational equations is selected. That order is (((C(N, M), S(N, M)), M=l,MMOD), 
N=l,NMOD), where coefficients may be omitted, but must not be out of order. SOR- 
DID therefore puts the IX array into the order expected by FIELD2 even if the user 
failed to do so himself. Following the sorting procedure, SORDID prints out the IX 
number and (for gravitational parameters) the location number in the input array as 
well as the values residing in the location when SORDID was called. SORDID is called 
before and after parameter estimation, so parameter values at these times are printed. 
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128 


144 ICO 176 


192 208 224 240 256 


Srp 

t , 1 

W 

GO 

*-■ 

S 9,l 

S 10,l 

s ll,l 

S 12, 1 

S 13, 1 

S 14.1 

S 15,l 

1 

? 7,2 

S 8,2 

S 9,2 

S 10,2 

S ll,2 

**12, 2 

S 13,2 

S 14,2 

S 15,2 

2 

S 7,3 

S 8,3 

S 9,3 

S 10,3 

S il,3 

S 1 2, 3 

S 13,3 

S 14, 3 

S 15, 3 

3 

S 7.4 

S 8,4 

6 9,.4 

S 10,4 

S ll,4 

S 12,4 

S 13,4 

S 14,4 

S 15,4 

4 

*7,5 

S 8,5 

S 9, 5 

S 10,5 

S ll, 5 

S 12,5 

S 13,5 

S 14,5 

S 15, 5 

5 

S 7,C 

S 8,6 

S 9, 6 

S 10,6 

S ll,6 

S 12,6 

S 13, 6 

• S 14, G 

S 15,6 

n 

U 

*7,7 









7 

C 8,3 

S 8,8 








8 


°9,9 

S 9,9 







9 



C 10,10 

8 10,10 






10 




C ll,ll 

S 11,H 





11 





C 12, 12 

S 12,12 




12 






C 13,13 

S 13,13 



13 







°14, 14 

S 14,14 


14 


C 15.15 S 15,15 


15 



SORDID FLOW CHART 
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SUBROUTINE SQRUT 


Calling Sequence: 

Purpose: 

Common Blocks Required: 
Subroutines Called: 

Input / Output 


CALL SQRUT (A,SA,N) 

SQRUT computes a square root of a positive 
semi-definite symmetric matrix. 

None 

None 


: : — — 1 

I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

w 

a' 

(N, N) 

Call List 

Input matrix, covariance type 

O 

SA 

N*(N+l)/2 

Call List 

Square root, packed upper- 
triangular 

i 

N 

1 

Call List 

\ ! 

1 

Matrix dimension 

i 


Method 

Given a matrix , A, we seek an upper-triangular matrix, S, such that 
T 

S S = A 


A recursive algorithm for computing S is 


s 

ii 



i=l,n 


s 

ji 


i < i 


i ( a - V' s s 
" ji E Jk ik 

k = l 


j=i + l, n 


The output matrix, SA, is "packed" as a singly-dimensioned array. 
SA = ( b iv s 12 , s 22 , s l3 , s 23 , s nn ) 
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SUBROUTINE TUMULT 


Calling Sequence: CALL TUMULT (A, B, C, NR, NC, IT) 

Purpose: TUMULT multiplies a packed upper- triangular 

matrix (or its transpose) times a rectangular 
matrix. 


Common Blocks Required: None 

Subroutines Called: None 

Input / Output 


r : 1 — — ~ \ 

I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

* 

Definition 

i 

A 

NR*(NR+l)/2 

Call List 

Packed upper-triangular matrix 

i 

B 

NR * NC 

Call List 

Rectangular matrix or vector 

0 

C 

NR * NC 

Call List 

Product matrix 

I 

NR 

1 

Call List 

Rows of B 

I 

NC 

1 

Call List 

Columns of B 

I 

IT 

1 

Call List 

_ n Tr ~ 1 0 : C= AB 

Transpose flag, IT = j 1 . c ^ A T B 


Discussion 

This subroutine uses the fact that A is upper-triangular to avoid unnecessary multiplications 
by zeroes in the product AB or A T B. A is packed as a singly-dimensioned array. 


A ( a ir a 12 , a 22 , a 13 , ..... 
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SUBROUTINE UPDATX 


Calling Sequence: 


CALL UPDATX 


Purpose: UPDATX adds the estimated state increment 

to the appropriate program arrays (estimated 
state) . 

Common Blocks Required: CNSDR, COMIC, DXCOM, FIELDM, IMCOM, 

INPUT, IXCOM, STATE 


Subroutines Called: None 


Input/ Output 


I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

I/O 

ELMI 

6 

COMIC (1) 

Estimated initial elements 

I/O 

DX 

N 

DXCOM(116) 

Incremental state 
(upper-loaded) 

I/O 

. 

X 

288 

FIELDM(l) 

Lunar gravitational field 
coefficients 

I 

ICNSDR 

20 

CNSDR (1) 

Considered parameter 
identifiers 

I/O 


1000 

INPUT(l) 

Input array 

m • 

- 

IX 

N 

IXCOM (1) 

State component identifiers 
(upper-loaded) 

I 

N 

1 

IXCOM(lOl) 

Number of incremental states 

I 

IGAIN , 

1 

IMCOM(7) 

Reciprocal gain factor 

I/O 

S(36) 

1 

STATE{36) 

Solar pressure coefficient 


68 











































Discussion 


Each incremental state component in DX is identified by the corresponding element 
of EX. The program arrays are incremented by DX with a unit gain as follows, 
accox*ding to IX. 

If IX(I) is DX(I) is DX is added to 

1^IX(I)*271 a gravitational coef- X [ IX( I ) ] in FIELDM 

ficient 

273 ^ IX(I) *278 an initial orbital element ELMI [lX(I)-272] in COMIC 

279 a solar pressure coef- X(36) in STATE 

ficient 

280 a solar pressure coef- A(198) in INPUT 

ficient 

A provision is made for adding the state increment to the estimate with an assigned 
gain (limit). If IGAIN is input as a number greater than unity, 

DX = DX/lGAIN 

Another provision permits exclusion of any incrementing of selected parameters. 

That is, if any IX(I) is found in the ICNSDR array of "considered' 1 parameters, then 
the state increment DX(I) corresponding to IX(I) will be set zero. This treatment 
has the effect of changing the weighting matrix without changing certain parameters 
during the estimation process. 

No flow chart is deemed necessary for this routine. 
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SUBROUTINE WLS 


Calling Sequence: 


CALL WLS 


Purpose: 


WLS is the driver for a weighted least 
squares estimation process. 


Common Blocks Required: CONST, DXCOM, IMCOM, INPUT, INTVRX, 

IXCOM, YCOM 

Subroutines Called: MODEL, PHOUSE, PUTINV, REDATA, SOLVE, 

SQRUT, TUMULT, UPDATX 


Input/Output 


I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

o 

QH 

5650 

DXCOM(201) 

Square root information matrix 

o 

HQY 

100 

DXCOM(lOl) 

Weighted residual vector 

o 

DX 

100 

DXCOM (1) 

Estimated state increment 

I/O 

ICVAR 

1 

INPUT(1096) 

Sensitivity computation type 
key 

I 

KMETH 

1 

INPUT(1036) 

Integration method key 

I 

NGROPT 

1 

INPUT(1064) 

Gradient treatment key 

I 

N 

1 

IXCOM(lOl) 

Number of state components 

I 

Ml 

1 

IXCOM(102) 1 Number of observations at 

f each measurement 

I 

M 

1 

IXCOM(103) 

Number of measurements 

0 

NQH 

1 

IXCOM(104) 

Number of elements in packed 
SRIM 

I 

IQ 

1 

IXCOM(105) 

. ... . „ 1. no weighting 

Weighting type flag 2. constant Q 
3. variable Q 

I 

IT MAX 

1 

IXCOM(106) 

Maximum number of iterations 
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I/O 

Symbolic 

Name 

Program 

Dimension 

Common 

Block 

Definition 

0 

LD 

1 

1 i IXCOM(107) 

Measurement number counter 

o 

IT 

1 

IXCOM(108) 

Iteration counter 

I 

' 

YHAT 

6 

YCOM(2) 

Estimated measurement counter 

I 

1 

IGAIN, LTOL 

1 

IMCOM(7, 8) 

Iteration convergence tolerance 
numbers 

I 

IM 

6 

IMCOM(l) 

Measurement selector array 

I 

Y 

6 

YCOM(9) 

Measurement vector 

I 

Q 

36 

YCOM(15) 

Measurement weighting vector 


Discussion 

WLS contains the gross logic for a weighted least squares estimation procedure, 
using a square-root information matrix formulation. The usual form for the WLS 
estimator is 



T -1 
H Q. 

1 1 



A 


dX 



where Y. is the measurement x'esidual vector (Y-YHAT) at the i-th measurement, 
-1 1 

Q. is the measurement weighting matrix, H. is the measurement sensitivity to 

1 a 1 

changes in the state and dX is the estimated state increment to be determined. In 

the square-root formulation the information matrix (bracketed term) is not formed 

per se. Instead, a square-root information matrix, W, is accumulated, where W 

obeys 

T V' T -1 

w w = L H i s H i • 

i=l 
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See the description of W ! s computation in subroutine PHOUSE. At each measurement, 


Q, 


1/2 


H is formed. This product is loaded into the QH array, beginning at location 
i 


NQH + 1 for use by PHOUSE in forming W. The weighted residual vector, HQY is 

accumulated as indicated by the following equation. 

m , -1/2 T , -l/2 y. \ 

HQY = ^ (Q. H.) (Q. Y.) 

i = l 


When HQY and W have been accumulated for M measurements, subroutine SOLVE is 
called to solve 

(W T W) dX = HQY 

for dX which then is used by UPDATX to increment the estimated state. 


WLS performs ITMAX iterations unless convergence, defined by 



i=l 


IGAIN 


10 


LTOL 


occurs first. During iterations after the first, the sensitivities {and, therefore, 

the square-root information matrix) may be held constant at their values on the 

_ 1 

NGROPT-th iteration. The weighted sensitivities, Q ? H, are written on logical 
23 as they are computed, then read back from 23 as they are needed. Of course, 
provision is made to avoid redundant computation of sensitivities by setting KVAR 
and NQ to zero. 


When less than 6 elements are used for estimation of state, the residual vector is 
upper-packed with the necessary Ml components and the weighting matrix is upper- 
left -packed as specified by the IM array. 
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WLS Flow Chart 


IT - 0 ( it = trial no . ) 



CALL TUMULT(SQ, DY, DY) Weight residuals, DY, 
CALL TUMULT(SQ, QH, QH> fand partials, QH 


Accumulate weighted 
residuals in HQY 

Accumulate square-root 
information matrix 


Compute estimated 

state increment 

Increment the estimated state 


CALL SOLVE \ 

i 

CALL UPDATX 




73 










m. 


THEORETICAL BACKGROUND 


3. 1 Definition of Symbols 


Capital letters will denote arrays (vectors or matrices) unless otherwise 
stated. Scalars will be denoted by lower case. A subscripted scalar will usually 
denote a particular element of an array. For example, W = [ w„ }. 


X, {x.} 
Y > 

Q* {%) 

n 

k 

m 

t 


E,{e.} 



d 


X 

A 

Y 

Y 

W '{ w ij} 

( ) 


State vector (gravitational coefficients, initial orbital elements, 
solar pressure coefficients) 

Measurement vector (orbital elements) 

Sensitivity matrix (of measurement to state) 

Measurement weighting matrix (error covariance matrix) 

Number of state elements to be solved for 

Number of observables in one measurement (usually 6) 

Number of measurements processed 

Time 

Orbital elements (measurement or a function of state and time) 
Zonal harmonic coefficients 

Tesseral harmonic coefficients in cosine of longitude 
Tesseral harmonic coefficients in sine of longitude 
An operator implying a variation of its operand 
Estimated state 

A 

Estimated measurement, Y(X,t) 

A 

Measurement residual, Y - Y 
Square root information matrix, upper triangular 

Superior numbers in parentheses refer to Bibliography 
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3. 2 Weighted Least Squares Parameter Estimation 

The fundamental weighted least squares estimation procedure will be derived 
with no attempt at rigor. 

Assume that variations, dX, in the state relate to deviations, dY, in the 
measurement by the linear relationship (3 . 1). 

dY = H dX (3.1) 

Presuming to have recorded the measurement vector, Y, and to have computed the 
measurement as a function of the state, Y(X), one may compute the measurement 
residual, Y. 

Y = Y - Y(X) ~ Y - Y (3.2) 

a 

If we were to hypothesize a change, dX, being made in the estimated state, we 
would expect that, under the linear assumption (3. 1), the residual would become : 

Y = Y - ( Y + H dX ) . (3.3) 

* ~T -1~ 

We seek that dX which minimizes the quadratic function, Y Q Y, where Q is a 

symmetric, positive -definite weighting matrix. One ready derivation of the 

A 

appropriate dX is found by differentiation of the quadratic function with respect to 
dX and equation of the result with zero, 

- 2 Y T Q _1 H = 0 (3.4) 

A. 

followed by solution of the resulting equation for dX. 

A T -1 -1 T -1 

dX = (HQ H) HQ (Y-Y) (3.5) 

T -1 

The matrix, HQ H, is called the "information matrix. " It must span the state 
space if it is to possess a unique inverse, although "pseudo-inverses" can always 
be computed to render solutions to (3.5). Iteration on (3.5) is performed by adding 

A A A 

dX to X and re-computing H and Y with the sum before re-applying (3.5), For 
certain problems, convergence can be obtained without re-computing H during the 
iteration. 
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3. 3 Square-Root Solution 


Certain advantages are to be obtained by using the "square root" form of (3. 5). 

The square root information matrix (SRIM) is defined here as the upper triangular 

matrix, W, such that (3.6) holds. 

W T W = H T Q _1 H (3-6) 

An upper triangular matrix is a square matrix characterized by having only zeroes 

T 

below its main diagonal. Conversely, a lower triangular matrix such as W has 
only zeroes above its main diagonal. The factorization of the information matrix 
can be accomplished by Cholesky decomposition. ^ In special cases, it can also 
be achieved by forming the' square root information matrix rather than the information 
matrix itself, so that the information matrix need never be factored. The solution 
of (3. 5) is obtained very efficiently, given W, by the method of Banachiewicz. 
Beginning with (3. 5) in the form of (3.7), 

(W T W) dX = H T Q~ 1 Y = D (3.7) 

we first solve for the dummy vector, B, of (3. 8) by forward substitution. 

W T B = D (3-8) 

A 

We then solve (3. 9) for dX by backward substitution. 


WdX = B 


(3.9) 


Denoting an element in the i-th row and j-th column of W by w.., 
substitution procedure is (3.10). 

\ = V w u 


me iorvvara 


b 2 = ,d 2 - W 12 b l )/W 22 


(3-10.1) 

(3.10.2) 


b. = (d. - 2- w..b. )/w.. 

3 3 £) U i ' 33 

The backward substitution procedure for n state variables is (3.11). 

dx = b /w 
n n nn 

dx = (b - w dx )/w , 

n-1 n-1 n-l,n n' n-l,n-l 


(3.10. i) 

(3.11.1) 

(3.11.2) 



(3.11.3) 
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This forward/backvvard substitution just described has the nice property that 

A. 

it provides dX satisfying (3.7) without explicitly inverting either the information 

matrix or its square root, thus saving the storage required for the SRIM and its 

inverse. If both the SRIM and its inverse are required for another reason or if 

storage is no problem, the substitution method is not quite as advantageous as 

simple inversion. If it is required, the inverse of the triangular SRIM is very simply 

( 2 ) 

obtained by another substitution procedure. ' For sequential information filtering, 
it is advantageous to retain the SRIM as opposed to replacing it in storage by its 
inverse. 
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3.4 Implementing the SRWLS Solution 


If successive measurement errors are assumed to be uncorrelated, the 

right-hand side of (3.7) can be accumulated as a simple vector sum and the SRIM 

( 1 , 2 ) 

can be accumulated by Householder reduction. ' 7 If we are trying to estimate 

n components of state and if we measure Ion quantities k-at-a-time, the dimensions 
of the component arrays of (3.7) are as shown in the following table. 

Table 3.1 Dimensions 
Array Dimension 


km x km 
km x n 


If successive measurements are uncorrelated by the assumed weighting matrix, then 


s 1 0 


• • Q_ 


where is a k x k matrix which weights components of the i-th measurement ag ains t 
each other and against Q. j ^ i. The Q. ^ are themselves positive definite and 
symmetric, possibly diagonal and possibly equal to each other. We can partition the 
residual into m k- vectors and H into m (k x n) matrices. If the residual vector for 
the i-th measurement is Y. and if H. is the k x n sensitivity matrix for the i-th 
measurement, the right-hand side of (3.7) is the simple vector sum (3.13). 


T -1 ~ w r r _-i ~ 

H Q Y =ZH 1 Q. 1 Y. 

fr, i i 1 

The information matrix can also be accumulated as a sum. 

T -1 V*\ ip _1 

HQ H = jT H. Q. H. 

£i 1 1 1 


(3.13) 
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The implemented square root method eliminates about half of the multiplications 
implied in (3.14) by accumulating the SRIM rather than the information matrix. If 
is the upper triangular SRIM after i-1 measurements have been included, the 
SRIM after the i-th measurement must obey (3.15). 


T T 

wrw. = w. ,w. . 

i i i-i i-i 


T -1 

+ H.Q. H. 

l l i 


Equation (3.15) can be written in partitioned form as follows. 


[w* p] 


W. 

i 

0 

L J 


T 

W i-1 


h t q : 1/2 

i i 


w 


i-1 


S 1/2h i. 


(3.15) 


(3.16) 


The computation of W. consists, then, of finding a means for converting 


w. , 


w. 

q- 1/2 h. 
1 1 


1 

0 


(3.17) 


where ^ and W. are n x n upper triangular matrices, where 


- 1/2 


H. and 0 are 
i 


kxn matrices and where W. satisfies (3.15). The process which has been implemented 

* - 1/2 - 1/2 
is called the Householder reduction algorithm. Since Q. and Q. H i must 

be formed for the HRA, the vector sum in(3.13) can be programmed as in equation (3.18). 


H 


V4 -.£( q: 1/2 h., t <^ 
/=/ 1 1 


- 1/2 


Y i> 


(3.18) 


The SRIM is initialized to zero to characterize the lack of a priori information. 
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Application to Lunar Field Determination 


The particular application of concern here is, of course, in determination 
of the harmonic coefficients of the lunar gravitational field. In this application, 
the state vector, X, consists of those coefficients we choose to solve for plus, 
optionally, the initial orbital elements (or Cartesian state) and one or more 
coefficients of the solar radiation pressure model. The state vector might have 
the following appearance. 


(example) 


X 





(3. 19) 


The measurement vector, Y, consists of elements of the lunar orbit. 

P 

e sin w 
Y _ e cos to 
M + to 
i 

a 


(3. 20) 


These elements are to be mean elements of the lunar orbit at specified times, referred 
to the Earth's mean equator and equinox of 1950. 0 and having been obtained from 
short-arc determinations using range and Doppler data. Sensitivities of Y to X are 
obtained by the Secant Method or, optionally, by integration of variational equations. 
Propagation of the estimated elements will be carried out by: 


1. choosing an initial estimate, X, of the state 

2. choosing initial elements, probably as the first measurement which would 
not then be processed as a measurement, 

3. integrating the averaged equations of motion across the time span of data. 
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Appendix A 

Derivation of Linear Variational Equations for 

the 

Method 7 Elements 

This appendix contains the expressions for the partial derivative matrices 
needed in the calculation of the rates for the linear variational equations for 
the method 7 elements. The method by which these equations are averaged 
is given at the end of the development. 


The equations of motion ( the planetary equations in Gauss's form) are 
E = F( E, C, t) = F4 E, r(E), C, t) 


or 


fv 


E = N{E) + M(E) 


R’ 

S’ 

W 1 


E = 


0 

q(-cos u) 
q sin u 
0 
0 

0 


"p 


0 

+ 

h 


0 


k 


0 


u 


vE|. 




r 


di_ 


0 


dt 




n 


0 







6x1 


6x1 


2qp (1 +k cos u + h sin u) 

r r r . 1 

’t 


(1 + — ) sin u + — h 
P P 


■ ) cos u + — k 
0 
0 

0 


M(E) 


6x3 


R 1 

S T 

W’ 


3x1 

0 


- q k — sin u cot i 
P 

r , 

q ~ h sm u cot l 
P 

r 

-q — sm u cot l 
H P 
r 

q — cos u 
P 

r sin u 


Q p sin i 


( 1 ) 


A _1 



We wish to describe the effects of small variations in the parameters C 
on the subsequent motion of the elements. Taking the variation of (1) with 


respect to the C's gives 

d 


f 

E 


SC 


d 

dt 


d E 

sc 


(/AjA + 14. 

-f£l 

a e 

I'-iE.l 

| \d E / dr 

d E ) 

d C 

\dc / 


where the normal brackets () represent explicit differentiation and the 
6-vector ^represents the cartesian position and velocity. The first 
column of the 6x6 matrix 


(#)- 


:ix (fr) lsgi 


smTs -7 


given by 
S N 


Sp W Sp 


. .. . , Sf , Sf 

and similarly for —7— and -rr- 

Sh /n ™\ Sk 


columns of the matrix 




The expressions for these first three 
are given in tables Al, A2, and A3. 


( 2 ) 


The disturbing acceleration vector, E' S’ W', is dependent upon the angular 
elements u, i, and D. as well as the time, viz: 


B’ 

S’ 


W’ 


[d] [q] [w] 

3x3 3x3 3x1 


where D is the transformation matrix from inertial coordinates to radial, 
circumferentialjOrbit. normal coordinates. The matrix D and its derivatives 
with respect to u, i, and 0 are given in table A7. The matrix Q represents 
the transformation from central planet body-fixed coordinates to the inertial 
integration frame. VV represents the disturbing acceleration due to the 
central planet non-sphericityand V is the familiar spherical harmonic representa- 
tion of the disturbing potential: 


00 


n 


V 


■2 £ 

n=l m=0 


n+1 


P ( sin ft) (C cos mX + S sxnmX) (3) 
n H nm nm 


where r, , and X are the body-fixed radius, latitude, and longitude of the 

spacecraft, B is the central planet's mean equatorial radius, and p m are the 

n 

associated Legendre polynomials. The gradient indicates differentiation with 


A - 2 



respect to central-planet-fixed coordinates. 


Because the matrix D is a function of i, Q, and u, the explicit derivatives in 
equation (2) must contain terms involving , -gy - ' , and ^ q . These 

expressions are given in tables A4, A5, and A6 respectively. 


The disturbing acceleration itself, VV, is treated as a function of the cartesian 
position vector and is differentiated in the usual way, viz: 


a W 
d r” 



T 


where r ' represents the cartesian position vector in a frame fixed in the 
central planet, r 1 is just the cartesian vector equivalent of the central body- 
fixed radius, r, latitude, /3 , and longitude A of the disturbing function given 
in equation (3). It is assumed that there are no velocity-dependent forces and 
the second term of (2) is 


mi°}\ 


5F \_ 

dr 


QJ /3 qT + J 3rd * 


SP 


where J is the Jacobian of the third body disturbing potential and J is 
3 i*d 

the gradient of the solar pressure disturbing force. Other disturbing forces 
should be treated in the same way. 


The remaining term in the curly brackets is the matrix representing the trans- 
formation from orbital elements to the cartesian position vector. Since 


, and r - l + h sin u + k cos u 


the first column of —77 is 

o B 




^ T 

l r = g>D 

r 

0 

T 

+ D 

~ar ~ 
dP 

BP 

0 

0 

. 0 _ 


and similarly for the other elements. The total matrix is given by 



A - 3 
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1 
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1 
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0 

1 

T 

— 
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0 


D 

0 

1 D 

0 

1 D 

0 

i 1 $ 1 

0 

a 

I 
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1 
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. 1 
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\ 
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0 

1 

1 

1 

0 

t 

1 

1 

0 

1 
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0 

i 
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0 
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0 

1 

0 

1 

0 

1 o ' 

0 

[ 

0 



0 

1 

0 

1 

0 

1 o 1 

0 

\ 
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‘ t 

. a d 

r 

T 

r r 2 I 

where X - 


+ D 

- — (h cos u - k sin u) 


0 


P 0 


0 


0 


The remaining term in equation (2) is the explicit derivative of the "plant" with 
respect to the parameters to be estimated. This differentiation is accomplished 
for the spherical harmonic coefficients by evaluation of the disturbing function 
with the desired coefficient set to unity and the other coefficients set to zero. 
This is possible because the harmonic coefficients enter the equations linearly. 
For efficiency purposes, this procedure is carried out internally in the sub- 
routine (FIELD2) that calculates the disturbing acceleration due to the central 
planet non- sphericity. The mechanics of this operation are explained in the 
FIELD2 subroutine description of Reference 2. 

Averaging the Linear Variational Equations 

In order to average the variational equations we must change variables in the 
formulation because the short-periodic parts of the derivatives of the fast 
variable with respect to h and k are not small. To change from u to H in the 
formulation of the linear variational equations, we simply calculate the two 
additional terms of the matrix M and their derivatives. 


In the first and second columns of the fourth row of M write 


M 


41 


- Vq ( k cos u + h sin u) - 2\Jl - e 2 ' ( ^ ) J 


A - 4 



and 



-# T 2 - 

V t‘ [e 


(k sin u - h cos u) (1 + — 


p>] 


in place of the zeros given earlier when the function 


was u instead of s, . 


The derivative matrices are modified similarly and we add to the tables A1 
through A 6 the following terms. 



- - - --- f Q-{ k cos u + h sin u) - 2 yl-e^ ( ~) 

2 L e p J 

pr [ — ( k sin u - h cos u) (1 + | )1 

\^p/ 42 2 ^ L e P J 


(tf ) 41 = VjF [l < f- >< k oos u + h sin u > + f sln u + ^pljTb 


- e‘ 


— ( - )1 
dh ^p'J 


( 


*m\ 

W 


42 




(-)(k sin u - h cos u)( 1 + — ) - cos u (1 + -) 
dh v e M p e P 


Q 


+ — (k sin u - h cos u) 

G 


± , 1,1 

dh k p 'J 


M\ _ Jp r <L , Q 

r k 


Q 

) ( k cos u + h sin u) + ^ cos u 


+ 2 ( F>/r^ 


- 


'] 


(tr) = #[£ ( ? ,lks,nu - hcbsu)< 

42 

“( £ )1 

dk v p ’] 


r . O 

1 + ~ ) + — sin u (1 + 
P e 


+ — { k sin u - h cos u) 
e 


(lf) 4l [; (hcosu ' ksin,1) ' 2 ^^ ^<p>] 


A - 5 


•73 |h 





where 


Q 1 _A . Q. ^ __d. ,Q, • _h 

e "^i_ e 2 +1 f dh e * de ' e ; 


VT ^2 dk ' e 


_d_ . Q. a d Q k 

de e yji _ e 2 


and , „ ^ ~ 2 

^(S) - - rV(?) 


de e 




2 'e 


Also 


-e‘ 

•* 

d , r- • . , r 

— (-) = - sinu ( -) 

dh v p ' p 


d . r . . r * 

-(-) = -COSU(-) 


dr r 

— (— ) = (k sin u - h cos u) ( — ) , 

du p p 


and 


” * f (P) so -“ { ^ ) * 0. 
p dp p 


Now, with the above modifications, we have written the linear variational eqviations 

for a system whose fast variable is 5 = <o+ M instead of u =n>+f but the right 

SS 

hand side is still in terms of u. The solution to these equations will involve 

>11 ^ 
rather than — — • and we must convert by writing 
S C 

q u _ *u_ _ _*H_ Wp 

5 C 3 C ?i H r^C nr^ 

whenever is needed. 


The only remaining modification is to change the term yjji. p/ r^ to yj P/ in 
accordance with the equations for the fast parts of the variables u and H respectively. 


The fourth equation now reads 


s = ^ 


s + V 

P 3 j=i 


M (RSW) 
4j J 


(4) 


A - 6 



and, if we say that n = 


Wa* , 


2 M. = ( V'' M (ESW) . I 

5C 5P ?C 5 h *c . dk bC 5C | 4 j Jj- 


The first three terms average to 

?n _ ?n ?£ + ^Jrf 5 h . 3~n p> k 
5C Vp a*h 5C ' ?)k $C 

and this large part of the variation of the fast variable can be added after the 
averaging integration with the help of 


^ n _ 3 _ n 

5 p " ' 2 (l-e*) 


A n 3 hn 5jn _ 3 kn 

?>h “ ~ ( 1 -e 2 ) , 5k " (l-e 2 ) . 


The variational equations then, are averaged in a way that is exactly analogous to 
the averaging of the method 8 variables. After the numerical averaging of the 
perturbation terms, the derivative of the average mean motion is added to the 
rate of the fourth variable and the averaging of the mean motion is taken outside 
the numerical averaging algorithm. 
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Table A 1 
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Table A3 




Table A 4 
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Table A 7 
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